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Abstract 

The  one-dimensional  motion  of  a  particle  is  analyzed  when  the  force  on 
it  is  inversely  proportional  to  its  displacement  and  directly  proportional 
to  the  elapsed  time.  Such  a  force  law  describes  a  projectile  in  a  musket 
barrel  that  is  propelled  by  a  hot  ideal  gas  where  either  the  number  of 
moles  or  the  temperature  increases  linearly  with  time  due  to  the  burning 
gunpowder.  A  particular  solution  to  Newton’s  second  law  is  found 
analytically  for  the  case  of  zero  initial  position  and  velocity.  For  more 
general  initial  conditions,  numerical  integration  is  used  to  find  the 
position  of  the  particle  as  a  function  of  time.  A  scaling  argument  shows 
that  at  long  times,  these  numerical  general  solutions  all  converge  to  the 
analytic  particular  solution.  Further  analysis  reveals  how  that 
convergence  occurs:  the  general  solutions  slowly  oscillate  about  the 
particular  solution  with  a  predictable  period  and  amplitude.  In  addition 
to  the  dynamics,  the  energetics  of  the  motion  are  analyzed. 


Basic  Dynamics  of  the  Ball  in  a  Musket 

A  lead  ball  of  mass  in  is  tamped  down  the  barrel  of  a  musket  of  cross- 
sectional  area  A,  so  that  it  rests  against  a  layer  of  black  powder  with  initial 
position  x0  =0  and  velocity  v0  =  0,  as  sketched  in  Figure  1.  Model  the 
system  by  making  three  assumptions  that  simplify  the  analysis  and  bring 
out  its  essential  features.  First,  suppose  [1]  the  powder  bums  at  a  constant 
rate  r,  creating  n  moles  of  hot  gas  at  temperature  T  as  a  function  of  time  t, 
so  that 


n  =  rt.  (1) 

Second,  assume  that  the  gas  expands  isothermally  [2]  as  the  ball  proceeds 
down  the  barrel,  so  that  T  is  constant.  Third,  neglect  the  losses:  Suppose 
that  there  is  no  friction  between  the  barrel  and  the  sliding  ball,  but  at  the 
same  time  assume  the  ball  fits  tightly  enough  that  there  is  no  leakage  of 
gas  around  it  [3].  (Historically,  lead  shot  would  often  be  wrapped  with 
linen  to  prevent  gas  from  escaping  while  minimizing  the  coefficient  of 


Summer  2012 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

2012 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2012  to  00-00-2012 

4.  TITLE  AND  SUBTITLE 

Analysis  of  the  Motion  of  a  Ball  in  the  Barrel  of  a  Musket 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

U.S.  Naval  Academy, Physics  Department , Annapolis, MD, 21402- 1363 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

Washington  Academy  of  Sciences  98.  Summer  2012 

14.  ABSTRACT 

The  one-dimensional  motion  of  a  particle  is  analyzed  when  the  force  on  it  is  inversely  proportional  to  its 
displacement  and  directly  proportional  to  the  elapsed  time.  Such  a  force  law  describes  a  projectile  in  a 
musket  barrel  that  is  propelled  by  a  hot  ideal  gas  where  either  the  number  of  moles  or  the  temperature 
increases  linearly  with  time  due  to  the  burning  gunpowder.  A  particular  solution  to  Newton?s  second  law  is 
found  analytically  for  the  case  of  zero  initial  position  and  velocity.  For  more  general  initial  conditions, 
numerical  integration  is  used  to  find  the  position  of  the  particle  as  a  function  of  time.  A  scaling  argument 
shows  that  at  long  times,  these  numerical  general  solutions  all  converge  to  the  analytic  particular  solution. 
Further  analysis  reveals  how  that  convergence  occurs:  the  general  solutions  slowly  oscillate  about  the 
particular  solution  with  a  predictable  period  and  amplitude.  In  addition  to  the  dynamics,  the  energetics  of 
the  motion  are  analyzed. 


15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

ABSTRACT 

18.  NUMBER 

OF  PAGES 

19a.  NAME  OF 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

17 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


36 


friction.)  Also  suppose  that  the  atmospheric  back  pressure  on  the  ball  can 
be  ignored  compared  to  the  forward  gas  pressure. 


x(t)  ) 

barrel  of  length  L 

G 

U  4 

powder 


'ball  of  mass  m 


Figure  1.  Ball  in  the  Barrel  of  a  Musket. 


The  gas  pressure  P  is  then  related  to  the  net  force  F  on  the  ball  by 
the  definition  of  pressure  as  force  per  unit  area,  P  =  F  /  A  =  ma 7  A  ,  where 
a  is  the  acceleration  of  the  ball  along  the  barrel  using  Newton’s  second 
law.  Furthermore,  the  volume  of  the  gas  is  given  by  V  =  Ax  where  x  is  the 
displacement  of  the  ball  (cf.  Figure  1),  and  thus  PV  =  max .  Treating  the 
gas  as  ideal,  i.e.,  PV  =  nRT  where  R  is  the  gas  constant,  it  follows  that 


xa  =  kt 

using  Equation  (1),  where  k  is  the  positive  constant 

,  rRT 

k  = - . 

m 


(2) 

(3) 


2  2 

Substituting  a  =  d  xt  dt  ,  Equation  (2)  becomes  a  nonlinear 
inhomogeneous  differential  equation  for  the  position  of  the  ball  as  a 
function  of  time,  x(t) .  A  particular  solution  of  it  can  be  immediately 
found  using  the  trial  form 

x  =  BtN  (4) 

where  B  and  N  are  constants  to  be  determined.  By  substituting  Equation 
(4)  into  (2)  and  equating  both  the  powers  and  prefactors  of  t  on  the  two 
sides  of  the  equation,  one  finds 

3  [4k 

N  =  -  and  B  =  J—,  (5) 

2  V  3 
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so  that 


x  =  (-|&)  t312 .  (6) 

The  first  and  second  time  derivatives  of  this  result  give  the  velocity  and 
acceleration  of  the  ball  as  a  function  of  time, 

v  =  kt  (7) 

and 


Equation  (6)  is  not  the  general  solution  of  Equation  (2)  because  it  does  not 
include  two  arbitrary  integration  constants.  Instead  it  is  the  particular 
solution  corresponding  to  x0  =  0  and  v{]  =  0  (which  happily  is  the  case 

most  relevant  to  a  ball  in  a  musket).  Later  in  this  article  we  will  consider 
how  to  find  solutions  for  other  initial  conditions  (see  the  section,  “General 
Solution  of  the  Differential  Equation,”  below). 

The  ball  is  in  the  barrel  during  the  time  interval  from  t  =  0  until 
some  later  time  t  =  tmax  .  Equations  (6)  and  (7)  then  imply  that  the  length 
of  the  barrel  is 

i  =  (9) 

and  the  muzzle  velocity  of  the  ball  is 

Hnax  —  \l^^t max.  •  (  ^  0) 

Equations  (9)  and  (10)  are  two  equations  in  two  unknowns  (k  and 
hnax)  if  L  anc*  Hnax  are  measured.  For  example,  the  58  Springfield  musket 
[4]  has  a  barrel  length  of  1  m  and  a  muzzle  velocity  of  290  m/s,  from 
which  one  deduces  that  k  =  5.4  km  /  s  and  /max  =5.2  ms .  Equations  (6) 

to  (8)  are  plotted  in  Figure  2  for  these  values.  Although  the  acceleration 
initially  diverges,  the  velocity  and  position  are  nevertheless  well  defined  at 
all  times.  The  rise  in  the  velocity  quickly  tapers  off,  demonstrating  that 
there  is  little  advantage  in  increasing  the  barrel  length  past  a  certain  point. 
(In  particular,  if  losses  were  included,  there  would  be  some  definite 
optimal  length  for  a  given  powder  charge.) 
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Figure  2.  Graphs  of  x  (in  m),  v  (in  km/s),  and  a  (in  m/ms2) 
versus  t  ranging  from  0  to  5.2  ms  if  k  =  5.4  km2  /s3  for  the  initial 
conditions  x0  =  0  and  v0  =  0  . 

Energetics  of  the  Powder,  Gas,  and  Projectile 

The  derivation  of  Equation  (2)  above  depends  on  the  fact  that  nRT 
increases  linearly  with  time  t.  Specifically,  it  was  assumed  that  T  is 
constant  and  that  n  =  rt .  Hot  gas  is  created  (starting  from  zero  moles)  by 
chemical  conversion  of  the  solid  powder.  For  simplicity,  consider  the  gas 
to  be  monatomic.  An  alternative  way  to  model  the  system  is  to  take  n  to  be 
constant,  while  T  =  rt,  i.e.,  the  temperature  is  no  longer  constant.  One 
could  now  think  of  the  gas  atoms  as  initially  existing  latently  in  the 
gunpowder  in  condensed  form  (which  classically  corresponds  to  absolute 
zero  temperature)  and  that  they  get  rapidly  warmed  up  by  thermal  energy 
transfer  Q  from  the  burning  charge. 

Reversible  thermodynamics  can  be  used  in  the  analysis  because  the 
gas  expansion  is  quasistatic  and  there  are  no  dissipative  losses.  Appendix 
A  shows  that  the  gas  is  always  in  quasi-equilibrium  by  treating  the 
expansion  of  the  gas  behind  the  ball  like  the  familiar  example  of  a  piston 
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in  a  frictionless  cylinder.  The  piston  slides  slowly  compared  to  the  speed 
of  sound,  provided  that  its  mass  is  much  larger  than  that  of  the  gas. 


n  =  constant  amount  of  gas  (initially  condensed) 


burning 

powder 


E  -  —nRT  internal  energy  of  gas 


Figure  3.  Relevant  energy  transfers  between  the  gunpowder, 
propelling  gas,  and  musket  ball  in  the  model  where  the  number 
of  moles  n  of  the  gas  is  taken  to  be  constant  and  its  temperature 
T  increases  linearly. 


The  energetics  of  the  musket  are  now  analyzed  by  reference  to 
Figure  3.  Consider  an  arbitrary  interval  of  time  between  0  and  t.  During 
this  interval,  the  monatomic  gas  ends  up  with  an  internal  energy  [5]  of 
E(t)  =  \.5nRT(t) ,  whereas  it  started  with  no  internal  energy  when 
condensed.  Thus  the  change  in  internal  energy  of  the  propellant  gas  is 


A E  = 


-nRT . 
2 


(11) 


During  this  time,  the  gas  does  work  W  on  the  ball,  calculated  from 

W  =  |  Fdx  =  m\  adx  (12) 

(or  equivalently  from  J  PdV .)  One  could  proceed  by  substituting 

Equations  (6)  and  (8)  to  convert  the  last  result  into  a  time  integral  that  can 
be  performed.  A  simpler  and  more  general  approach  is  to  use  the 
definitions  of  velocity  and  acceleration  to  rewrite  Equation  (12)  as 

W  =  m^vdt  =  ^mA([V2))  =  AK  (13) 

where  K  =  \mv~  is  the  kinetic  energy  of  the  ball.  Equation  (13)  is 
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simply  a  derivation  of  the  work-kinetic-energy  theorem  of  the  ball  since  it 
is  treated  here  as  a  particle.  For  the  particular  solution  of  interest  (which  is 
also  the  asymptotic  solution  for  large  times  for  any  choice  of  initial 
conditions,  as  is  proven  in  the  section  below,  “Scaling  Argument  to  Find 
the  Asymptotic  Behavior”),  Equations  (1),  (3),  and  (7)  can  be  substituted 
into  this  result  to  obtain 

W  =  —nRT .  (14) 

2 

Finally,  the  first  law  of  thennodynamics  applied  to  the  gas  implies  that 

A E  =  Q-W  =>  Q  =  3nRT  (15) 

using  Equations  (11)  and  (14).  In  words,  these  results  prove  that  energy  is 
being  transferred  from  the  powder  to  the  gas  and  projectile  at  the  constant 
rate  dQ  /  dt  =  3 nRr  .  Half  of  that  added  energy  goes  into  warming  up  the 
gas,  increasing  its  internal  energy,  while  the  other  half  goes  into 
accelerating  the  ball,  increasing  its  translational  kinetic  energy.  Thus  the 
energy  transfer  efficiency  from  the  powder  to  the  ball  is  50%  even  for  this 
ideal  musket. 


General  Solution  of  the  Differential  Equation 

In  this  section  we  explore  the  nonlinear  differential  equation 
xa  =  kt  from  a  more  general  mathematical  point  of  view.  The  physical 
application  presented  above  helps  interpret  these  general  numerical 
results. 


Reparameterizing  the  Equation 

Equation  (2)  is  second  order  and  so  its  general  solution  must 
contain  two  constants  of  integration,  corresponding  to  the  initial  position 
A'o  and  initial  velocity  v0,  in  addition  to  the  adjustable  value  of  k.  In  total, 
there  are  thus  three  parameters  in  the  family  of  solutions.  However,  some 
of  these  parameters  can  be  eliminated  by  rescaling  the  units  of  length  and 
time.  This  reparameterization  is  performed  in  different  ways,  depending 
on  the  initial  conditions. 

Choose  the  +x  axis  to  point  in  the  direction  of  motion  of  the 
particle  at  long  times.  Then  for  k  >  0 ,  the  position  of  the  particle  can 
never  be  negative.  Appendix  B  discusses  what  happens  if  the  ball  is 
launched  toward  the  origin  and  shows  that  x  cannot  reach  (or  cross)  zero 
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no  matter  how  large  the  initial  speed  is.  The  particle  is  repelled  from  the 
origin  and  must  eventually  move  in  the  positive  x  direction,  either  because 
it  initially  started  moving  in  that  direction  or  because  it  reversed  direction 
after  exactly  one  bounce. 

Therefore,  at  times  t  other  than  zero,  the  position  x  must  always  be 
strictly  positive.  In  contrast,  at  t  =  0  the  initial  position  x{)  can  be  either 
positive  or  zero.  For  example,  Equation  (2)  indicates  that  x0  can  be 
nonzero  if  the  initial  acceleration  is  zero,  whereas  Equation  (6)  applies  to 
the  special  case  x0  =  0  .  If  x0  =  0 ,  then  the  initial  velocity  v0  cannot  be 

negative,  because  otherwise  the  particle  would  thereafter  move  in  the 
direction  of  negative  x,  contrary  to  our  choice  of  axis.  The  various  possible 
initial  conditions  can  therefore  be  divided  into  three  classes. 

Class  A:  xn  is  positive  and  vn  is  arbitrary 
Define  a  characteristic  length  x  =  x0  and  a  characteristic  time 
t  =Xq  k  .  These  two  quantities  can  be  used  to  define  a  characteristic 

speed  v  =  x/t=x0  k  and  a  characteristic  acceleration 

/~2  —1/3  2/3 

t  =Xq  k  .  In  effect,  the  units  of  distance  and  time  have  been 
chosen  to  scale  away  x0  and  k,  reducing  the  three-parameter  family  of 
solutions  to  a  form  that  depends  only  on  v0. 

Class  B:  Xn  is  zero  but  Vr,  is  positive 

3  —1 

Now  define  the  characteristic  length  as  x  =  v{]  k  and  the 
characteristic  time  as  t  =v^k  .  Then  the  characteristic  speed  is 
v  =  x  / 1  =v0  and  the  characteristic  acceleration  is  a  =  x/t2=  v^k  .  That 
is,  the  units  of  distance  and  time  have  been  chosen  to  eliminate  l)0  and  k. 
Given  that  x0  =  0  ,  the  family  of  solutions  to  Equation  (2)  reduces  to  zero 
adjustable  parameters  with  two  specified  initial  conditions. 

Class  C:  both  x0  and  ih  are  zero 

In  this  case,  there  is  only  one  parameter  in  the  original  problem 
and  so  we  cannot  independently  define  both  a  characteristic  length  and 

time.  Arbitrarily  choosing  x  =  1  m  in  base  SI  units,  then  t  =  x“  k  , 

~  1/3  1/3  —1/3  2/3 

v  =  x  k  ,  and  a  =  x  k  .  This  class  corresponds  to  the  particular 
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solution  already  found  in  the  first  section  of  this  paper,  “Basic  Dynamics 
of  the  Ball  in  a  Musket.” 

For  any  of  the  three  classes,  dimensionless  kinematic  variables  can 
be  introduced  as  t'  =  t  It ,  x  =  x  /  x ,  v  =v  /  v ,  and  a'  =  a/  a.  In  terms  of 
these  dimensionless  variables,  Equation  (2)  becomes 

x'a'  =  t'  (16) 

where  v  =  dx  l  dt'  and  a  =dv'  /  dt' .  The  initial  conditions  for  class  A  are 
x'0  =  1  and  Vq  =  Xq ll3k~ll3v0  .  The  differential  equation  has  thus  been 
recast  in  a  form  that  only  depends  on  the  single  combined  value  v'0 .  That 
makes  it  easier  to  investigate  and  plot  its  family  of  solutions.  For  class  B, 
the  initial  conditions  are  uniquely  specified  as  x'0  =  0  and  v'q  =  1 .  Using 
l’Hopital’s  rule,  the  initial  acceleration  is  then  a'Q  =  1 . 

Numerical  Solution 

Equation  (16)  has  no  discernible  closed-form  analytic  solution  in 
general,  but  one  can  numerically  integrate  it.  Different  methods  can  be 
used  for  this  purpose,  depending  on  the  desired  accuracy  and  ease  of 
calculations.  To  start  with,  the  first-order  Euler-Cromer  method  [6]  can  be 
implemented  in  a  spreadsheet.  Given  values  x'  and  v'  at  any  time  t' , 
their  values  at  the  next  time  step  t'M  =  t\  +  At'  (where  At' =  0.1  say, 
corresponding  to  a  time  step  in  real  units  of  (7 10)  are  sequentially 
calculated  as 

v'+\  =»'  +  a'At'  =  v'  +  (t'/  x' )  At' 

x'i+i  ~  x'i  +  v'i+iAt' 

with  starting  values  for  class  A  of  t'Q  =  0 ,  x'()  =  1 ,  and  any  chosen  value  of 
v'q  .  The  results  are  plotted  as  the  solid  curves  in  Figure  4  for  four  values 
of  v'q  .  A  similar  numerical  integration  for  class  B  (when  x'0  =  0 ,  v'q  =  1 , 
and  a'0  =  1 )  results  in  a  curve  that  is  almost  indistinguishable  from  the 
dashed  curve  (corresponding  to  class  C). 

The  equation  for  xi+l  involves  the  updated  velocity  v'+]  rather  than 
the  previous  value  v\ ,  in  contrast  to  the  equation  for  v'+]  which  uses  the 
previous  value  of  the  acceleration  a\ .  That  is  the  hallmark  of  Cromer’s 
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modification  of  the  standard  Euler  method,  classified  as  a  symplectic 
integrator.  Symplectic  methods  are  the  preferred  choice  when  the 
Lagrangian  L  has  no  explicit  time  dependence,  but  they  can  be  used  even 
when  L  is  time  dependent  (as  is  shown  in  Appendix  C  to  be  the  case  here). 
By  comparison  with  the  results  of  the  second-order  leapfrog  integration  in 
Appendix  D,  the  Euler-Cromer  method  obtains  values  of  position  and 
velocity  that  are  found  to  be  accurate  to  0.4%  or  better,  so  Equation  (17) 
suffices  to  generate  Figure  4. 


dimensionless  time 


Figure  4.  Graphs  of  x  versus  t' .  The  four  solid  curves  were 
calculated  numerically  using  Equation  (17)  for  x'()  =  1  and  v'0 
equal  to  2,  1,  0,  and  -1  (from  top  to  bottom  at  t'  =  1 ).  The  dashed 
curve  is  a  plot  of  Equation  (18). 


Scaling  Argument  to  Find  the  Asymptotic  Behavior 

The  dashed  curve  in  Figure  4  is  a  plot  of  the  solution  represented 
by  Equation  (6)  in  dimensionless  form, 


4  / 3/2 
3 


(18) 


which  is  the  particular  solution  of  Equation  (16)  for  class  C,  i.e.,  initial 
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values  x'q  =  0  and  v'0  =  0  .  To  investigate  how  Equation  (18)  relates  to  the 
numerical  results  in  the  “Numerical  Solution”  section  above,  introduce 

new  variables  x"  =  x  /  s3  and  t"  =  t'/  s2  where  s  is  an  arbitrary  positive 
scaling  factor  that  zooms  Figure  4  in  or  out  (albeit  not  equally  for  the  two 

axes).  Defining  if  =  dx"  /  dt"  and  a"  =  d2x"  /  dt"2 ,  Equation  (16)  becomes 

x"a=f  (19) 

with  initial  conditions  Xq  =  Xq  /  s3  and  Vq  =v'Q  Is.  Then  if  5  — >  00 ,  the 
solution  corresponding  to  Xq  =  0  and  Vq  -  0  is  obtained,  namely 

x'  =  (4/3)1/2fi'3/2.  Therefore  Equation  (18)  is  an  asymptotic  solution  of 
Equation  (16)  at  large  times  for  any  initial  values.  To  verify  that 
conclusion,  Figure  4  was  replotted  with  the  range  of  t'  values  increased 
25-fold.  On  this  expanded  scale,  the  solid  curves  are  indistinguishable  by 
eye  from  the  dashed  curve. 

Oscillations  about  the  Asymptotic  Solution 

To  quantify  the  manner  in  which  the  general  solutions  approach 
the  asymptotic  curve,  dimensionless  residuals  Ax'  are  computed  by 
subtracting  the  analytic  particular  solution  given  by  Equation  (18)  from 
any  solution  x'(t')  computed  numerically.  Whenever  a  numerical  solution 
has  a  smaller  value  of  x  (at  a  given  t' )  than  the  asymptotic  solution,  it 
has  a  larger  acceleration  ( i.e .,  second  derivative),  and  vice-versa, 
according  to  Equation  (16).  The  numerical  solution  thus  repeatedly 
catches  up  to  and  crosses  the  asymptotic  curve,  oscillating  about  it.  It  is 
found  that  these  oscillations  have  both  increasing  period  and  increasing 
absolute  amplitude. 

A  physical  explanation  for  these  oscillations  can  be  found  in  the 
ballistic  situation  of  the  first  section,  “Basic  Dynamics  of  the  Ball  in  a 
Musket.”  If  the  ball  gets  ahead  of  the  position  it  has  in  the  particular 
solution  at  the  same  time,  the  gas  becomes  under-pressured  in  comparison 
and  the  ball’s  acceleration  drops.  That  pennits  the  gas  pressure  to  “catch 
up,”  but  the  inertia  of  the  ball  leads  to  an  overcorrection,  and  the  cycle 
now  reverses. 

Returning  to  the  residuals,  the  amplitude  of  Ax'  is  proportional  to 
as  t'  — >  00  ,  according  to  the  analysis  in  Appendix  D.  Therefore  the 
relative  amplitude  of  the  oscillations,  Ax'/x' ,  decreases  to  zero  in 
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proportion  to  1  It'  for  large  times.  In  this  relative  sense,  solutions  of 
Equation  (16)  for  any  initial  conditions  converge  onto  the  asymptotic 
solution,  Equation  (18). 


Conclusions 

The  force  law  described  by  xa  =  kt  has  a  simple  form  but  exhibits 
rich  behavior.  Physically,  it  describes  the  idealized  motion  of  a  ball  in  an 
unrifled  musket.  Another  situation  that  gives  rise  to  Equation  (2)  is  the 
radial  motion  of  a  particle  (of  mass  m  and  charge  q)  in  the  electric  field  of 
a  long  straight  wire  whose  linear  charge  density  is  proportional  to  time, 
A-at .  Then  k  =  aq  /  2 ne0m  (where  e0  is  the  pennittivity  of  free  space), 

which  allows  the  possibility  of  a  negative  value  of  k  that  the  ballistic 
application  does  not.  Investigation  of  the  motion  for  k  <  0  could  make  an 
interesting  student  project.  Further  study  of  the  equation  a^ttx  may 
uncover  additional  applications  and  intriguing  behavior. 

Mathematically,  the  solution  of  Equation  (2)  involves  power-law 
behavior  (of  the  particular  solution),  oscillatory  behavior  (of  the 
residuals),  and  exponential  behavior  (of  the  intervals  between  zero 
crossings  of  the  oscillations).  It  calls  for  a  diverse  combination  of  analysis 
techniques,  including  insights  from  physics,  trial  solutions  of  differential 
equations,  scaling  laws,  graphical  methods,  algebraic  approximations,  and 
computational  algorithms. 


Appendix  A:  Quasistatic  Expansion  of  the  Gas 
Pushing  the  Ball  in  the  Barrel 


At  long  times  (if  not  initially),  the  speed  of  the  musket  ball  is  given 
by  the  slope  of  the  dashed  curve  in  Figure  4,  obtained  by  substituting 
Equation  (3)  into  (7)  to  get 


l3nRT 
V  mball 


(A.l) 


after  replacing  n  =  rt  from  Equation  (1).  Here  the  subscript  “ball”  has 
been  added  to  v  and  m  to  emphasize  that  those  two  variables  refer  to  the 
speed  and  mass  of  the  projectile  specifically.  On  the  other  hand,  the  root- 
mean-square  speed  of  the  atoms  in  a  monatomic  ideal  gas  is  known  from 
kinetic  theory  or  equipartition  [5]  to  be 
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V, 


3  nRT 


gas 


m 


gas 


(A.2) 


where  mgas  is  the  total  mass  of  the  gas  in  the  musket.  (The  speed  of  sound 
has  almost  the  same  value,  obtained  by  replacing  the  factor  of  3  in  this 
equation  by  the  adiabatic  exponent  y=5/3.)  Comparing  Eqs.  (A.l)  and 
(A.2),  one  sees  that  the  speed  of  the  ball  will  always  be  small  compared  to 
typical  molecular  speeds  provided  that  the  mass  of  the  ball  is  much  larger 
than  the  total  mass  of  the  gas.  In  that  case,  the  gas  expansion  is  said  to  be 
quasistatically  slow. 


Appendix  B:  Motion  of  the  Ball  for  Negative  Initial  Velocities 

One  cannot  compress  the  volume  of  the  gas  in  the  barrel  to  zero  (at 
nonzero  temperatures).  However,  if  one  runs  the  Euler-Cromer  simulation 
specified  by  Equation  (17)  with  say  v'Q  =-3  and  At'  =  0.1 ,  the  trajectory 

appears  to  cross  the  horizontal  axis  (near  t'  =  0.34 )  and  become 
increasingly  negative  thereafter.  This  zero  crossing  is  an  artifact  of 
inaccurate  numerical  integration.  It  only  occurs  when  the  time  step  is  large 
enough  that  the  simulation  can  “jump”  over  the  divergence  in  the 
repulsion  that  occurs  at  x  -  0 .  We  can  prove  that  such  a  jump  does  not 
occur  for  infinitesimally  fine  time  steps  as  follows. 

Consider  the  specific  work,  i.e.,  work  per  unit  mass,  |  a'dx  as  the 

ball  approaches  x  =  0 .  Then  t'  is  approximately  constant  (for  example  at 
about  0.34  if  v'0  =-3)  during  the  small  time  interval  that  x'  is  smaller 
than  some  initial  value  x' ,  say  0.1.  Now  according  to  the  work-kinetic- 
energy  theorem, 


I)'2  -  V: 


(2  =  2  [  a'dx'  ~  2t'  In  ^ 

d  V. 


(B.l) 


using  Equation  (16).  The  right-hand  side  of  this  equation  is  negative 
(because  x'  <  x[  )  and  thus  the  ball  slows  down  as  it  approaches  x  -  0 .  But 
it  can  never  reach  x  -  0  because  the  left-hand  side  can  never  get  smaller 
than  -v'2  in  value.  Consequently,  even  if  the  negative  initial  velocity  is 

very  large  in  magnitude,  the  ball  will  necessarily  bounce  off  the  (infinite) 
potential  barrier  at  the  origin,  as  occurs  for  the  curve  corresponding  to 
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v'Q  =  -1  in  Figure  4. 


Appendix  C:  Canonical  Mechanics  of  the  System 

An  explicitly  time-dependent  Lagrangian  L  and  Hamiltonian  H  can 
be  constructed  as  follows.  Equation  (B.l)  suggests  a  potential  energy  U(t) 
that  is  logarithmic  in  the  position  and  thus 

L  =  K-U  =  j  mx2  +  mkt  In  x  (C .  1 ) 

where  x  =  v.  The  momentum  conjugate  to  the  position  x  is 
dL  . 

p  =  —  =  mx  (C.2) 

dx 

and  thus  the  Hamiltonian  is 


1  2 

H  =  px-L  =  \mx  -  mkt  In  x . 


(C.3) 


Then  the  equation  of  motion  is  obtained  from  Hamilton’s  equation  as 


dH  ..  mkt 

p  = - =>  m  x  = - 

dx  x 


(C.4) 


which  rearranges  into  xa  =  kt .  Alternatively,  this  equation  of  motion  can 
be  obtained  from  Equation  (C.l)  using  the  Lagrange  equation 


d_dL_dL_ 
dt  dx  dx 


(C.5) 


Appendix  D:  Harmonic  Oscillations  of  the  Residuals 

Dropping  the  primes  so  as  to  unclutter  the  notation,  the 
dimensionless  force  law  is 


..  t 
x  =  — 

x 

from  Equation  (16),  with  a  particular  solution  of 


x. 


=  P3'2 


(D.l) 


(D.2) 


Define  the  residual  Ax  as  the  difference  between  a  general  solution  of 
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Equation  (D.l)  for  x  and  the  right-hand  side  of  Equation  (D.2). 
Examination  of  numerical  results  from  Equation  (17)  suggests  that  the 
residual  (for  any  value  of  the  dimensionless  initial  velocity  v0 )  oscillates 
harmonically  in  the  logarithm  of  the  elapsed  time  t  with  an  amplitude  that 
increases  in  proportion  to  the  square  root  of  the  time.  To  verify  this 
suggestion,  let  the  scaled  residual  z  be  defined  as  Ax  divided  by  yft ,  so 
that 


z  =  xt 


-1/2 


which  can  be  solved  for  x  to  get 
x  =  zt1/2+J±t3/2. 


(D.3) 


(P-4) 


Take  the  second  derivative  of  this  equation  with  respect  to  time,  and 
substitute  both  it  and  Equation  (D.4)  into  (D.l)  to  obtain 

[l+^liz/t)  =  l+^(zt+z~iz/t)-  (P-P 


Since  it  will  be  shown  that  the  amplitude  of  the  scaled  residual  z  levels  off 
in  value  asymptotically  (cf.  Figure  5),  z/t  must  approach  zero  for  large  t. 
Thus  the  left-hand  side  of  Equation  (D.5)  can  be  approximated  using  the 
binomial  expansion  to  first  order.  The  result  can  be  rearranged  to  get 

zt2  +  zt  =  -\z.  (D.6) 

The  left-hand  side  of  this  equation  can  be  identified  as  the  logarithmic 
second  derivative 

d2z  t 

- T  =  zt2+zt .  (D.7) 

(dint)2 

Therefore  Equation  (D.6)  implies  that  the  (logarithmic)  second  derivative 
of  z  is  proportional  to  -z,  characteristic  of  simple  harmonic  motion. 
Consequently  the  residual  oscillates  semi-logarithmically  with  a  period  of 
2x^2  ,  /. e. ,  a  half-period  corresponds  to  a  ratio  of  adjacent  zero-crossings 
of  z(t)  that  equals  exp^W^  ~  85.02  . 
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Figure  5.  Semi-logarithmic  plots  of  the  scaled  residual  Ax'/-J7 
against  the  dimensionless  time  t'  for  the  three  indicated  values 
of  v'0  with  Xq  =  1  . 


Figure  5  plots  the  scaled  residual  z  as  a  function  of  time  on  semi- 
logarithmic  axes  for  x0  =  1 .  The  three  curves  correspond  to  different 

values  of  v0  between  0  and  -1.  Zero  crossings  for  these  three  curves  are 
listed  in  Table  1.  For  large  t,  the  ratio  between  successive  zero  crossings  is 
in  excellent  agreement  with  the  asymptotic  value  85.02  predicted  by 
Equation  (D.6). 

The  results  in  Table  1  require  x  values  accurate  to  better  than  1  part 
in  109,  because  z  is  a  small  difference  between  large  numbers.  To  achieve 
that  level  of  accuracy,  a  C++  program  [7]  was  written  to  calculate  the 
residuals  over  a  longer  range  of  times  and  with  higher  accuracy  than  can 
be  done  using  the  spreadsheet  solution.  The  program  replaces  the  first- 
order  calculations  of  Equation  (17)  with  the  second-order  leapfrog 
calculations 

xi+\  =xi+vAt  +  -  —  Of  and  vM=vi  +  ~  ^  +  —  At.  (D.8) 

2Vxi)  2Vxi  xi+ 1) 
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An  adaptive  time  step  size  is  used:  At  is  gradually  increased  as  t  increases, 
to  keep  pace  with  the  exponential  increase  in  the  period  of  oscillations.  An 
even  higher-order  symplectic  integrator  [8]  was  used  as  a  further  check  on 
the  results. 


Table  1.  Zero-crossing  times  in  Figure  5  and  their  ratios,  accurate  to 
1  part  in  1 0^,  for  three  different  values  of  the  dimensionless  velocity. 
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